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University of Southern California 

In this paper we lay the foundation for a numerical algorithm 
to simulate high-dimensional coupled FBSDEs under weak coupling 
or monotonicity conditions. In particular, we prove convergence of a 
time discretization and a Markovian iteration. The iteration differs 
from standard Picard iterations for FBSDEs in that the dimension of 
the underlying Markovian process does not increase with the number 
of iterations. This feature seems to be indispensable for an efficient 
iterative scheme from a numerical point of view. We finally suggest 
a fully explicit numerical algorithm and present some numerical ex- 
amples with up to 10-dimensional state space. 

1. Introduction. Motivated by the aim to simulate iiigh-dimensional cou- 
pled forward-backward stochastic differential equations (FBSDEs) we study 
a time discretization and a Markovian iteration for equations of the form 



(1.1) 



Xt = x+ f\{s,Xs,Ys)ds+ f a{s,Xs,Ys)dW, 

Yt = g{XT)+ I f{s,Xs,Ys,Zs)ds- f Z^dW 
Jt Jt 



where 6, cj, /, g are deterministic and Lipschitz continuous functions of lin- 
ear growth which are additionally supposed to satisfy some weak coupling 
or monotonicity condition. The solution consists of a triplet (X, y, Z) of 
adapted processes, which are called the "forward part," the "backward part" 
and the "control part" respectively. The presence of the control part Z is 
crucial to find a nonanticipative solution. It often has an intuitive interpre- 
tation, for example, as an investment strategy in financial applications; see 
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El Karoui, Peng and Quenez [12]. Note that (1.1) is not in its most general 
form, since Z does not couple into the forward SDE. 

Most of the numerical algorithms for coupled FBSDEs, with the notable 
exception of Delarue and Menozzi [10], exploit the relation to quasi-linear 
parabolic PDEs via the Ma-Protter-Yong four-step-scheme [17]. Under ap- 
propriate conditions, (X, y, Z) are connected by 



(1.2) Yt = u(t,Xt), Zt = v{t,Xt)^Ux{t,Xt)a{t,Xt,u{t,Xt)), 



The main focus in these approaches is on the numerical solution of the 
PDE (1.3); see Douglas, Ma and Protter [11], Milstein and Tretyakov [20] 
and Ma, Shen and Zhao [18]. Since the PDE approach requires existence 
of a classical solution to (1.3), there is typically need for some smoothness, 
boundedness and regularity conditions, such as uniform ellipticity of the 
differential operator. For low-dimensional problems, under such regularity 
conditions, the PDE approach may generally be regarded as superior to 
Monte Carlo simulation concerning accuracy and speed. However, solving 
(1.3) numerically by standard PDE techniques becomes more difficult, if 
not impossible, with increasing spatial dimension. Hence, it seems necessary 
to tackle the FBSDE (1.1) directly by probabilistic means in order to solve 
(1.1) numerically in situations which are beyond the limitations of the PDE 
approach. 

A natural time discretization of equation (1.1) is 



where h = ^ and ti = ih,i = 0,1, ... ,n, and AP1/j+i = VFt-^-^ — W^. . Here, 
of course, E^. denotes the conditional expectation E{-\J^t.}. This time dis- 
cretization was investigated in detail by Zhang [25] for decoupled FBSDEs. 
Note that Z in (1.2) and in (1.4) may be considered analogous to each 
other. Indeed, the expression for Z can be rewritten as the Malliavin deriva- 
tive DtYt of Y and, applying integration by parts under the conditional 
expectation, 




(1.4) 




[ Y^ ^ Eu{Yr+i + f{k,Xf,Yr+r,Zm, 
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is a natural discretization of the Malliavin derivative of Y^. Concerning 
the iteration (1.4), it is crucial to notice that X is discretized forwardly 
and Y is discretized backwardly. Hence, (1.4) is by no means an explicit 
discretization in the present situation due to the coupling and therefore 
cannot be implemented directly. Note, however, that one can rewrite 

(1-5) y,"=<(xr), z^=vnxr), 

where 



(1.6) 



X^_lY = x + biti,x,u^{x))h + a{t^,x,uf{x))AW^+l, 
I <(x) ^ EiYl"^^ + /(t„x,y,:;f ,<(x))M. 



Equation (1.6) is still implicit in u", but truly backward in time. Com- 
bined with a local updating technique, it serves as starting point for the 
probabilistic scheme in Delarue and Menozzi [10]. This type of scheme re- 
quires, however, apart from estimating the expectations, a discretization of 
the state space. Such space discretization may again become prohibitive, 
when the dimension increases. 

We, hence, propose to combine the time discretization (1.4) with an it- 
erative scheme. It is known from results by Antonelli [1] and Pardoux and 
Tang [22] that, under weak coupling or monotonicity conditions, (1.1) has a 
unique solution {X, y, Z) which can be constructed via a Picard iteration 

Xr = x+ f\{s,XT,Yr-')ds+ f'a{s,XT,Yr'-')dWs, 
Jo Jo 

Yr = gim + [ f{s, K\ Yr, zd ds - z^ dw^, 



(1.7) 



starting at Y^ = 0. The drawback of (1.7) is that the dimension of the 
underlying Markovian process increases with the number of iterations m. 
Precisely, one can easily see that Y^ is a function of X^ and, hence, the 
right-hand side of the SDE for X^ depends on X^ (through Y^) and X"^ . 
Proceeding this way, one observes that generally is not Markovian, but 
only the extended system {X^ ■, ■ ■ ■ ,X'") is. Consequently, Yr is a function 
of time and (X^, . . . , X*"), and therefore, the computational effort to 
estimate rapidly increases with the number of Picard iterations. This 
renders a combination of (1.4) with a Picard iteration like (1.7), which was 
recently suggested by Riviere [23] in theory, impractical from a numerical 
point of view. The stochastic control approach in Cvitanic and Zhang [8], 
which iterates over Z, faces the same kind of difficulty. 
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In this paper we introduce an alternative iteration in a way that the di- 
mension of the underlying Markovian process does not change in the number 
of iterations. It reads, in discretized form, = 0, and 

— X, 

Xi^^=Xi' +b{ti,X^' ,n.' (X,' ))h 

I /j. v'n,m n,m— 1 / x^n,m\\ A TiA 



\^n,m A „/' vn,m\ 



The main advantage is that here yj"'"^ is a function of time and X"'"*, but 
does not depend on (X"''^, ;U = l,...,m — 1). Establishing the convergence of 
this new "Markovian" iteration turns out to be more involved than for the 
standard Picard iteration, because controlling the Lipschitz constant and 
the linear growth of ^"'"^(x) uniformly in i,n,m becomes crucial. This is 
indeed the reason why we cannot allow Z to couple in the forward SDE at 
the current state of our research. 

We also indicate how this discretized Markovian iteration may be trans- 
formed into a viable numerical scheme, replacing the conditional expec- 
tations by simulation based least squares regression and estimating u"'"* 
this way. Such an estimator was introduced by Carriere [6], Longstaff and 
Schwartz [16] and Clement, Lamberton and Protter [7] in the context of 
American options and is applied by Gobet, Lemor and Warin [14] and Ben- 
der and Denk [3] for decoupled FBSDEs. Although a convergence analysis for 
this estimator in the present context of a coupled FBSDE is beyond the scope 
of this paper, we illustrate by some examples with up to 10-dimensional state 
space that the proposed numerical algorithm works in practice. 

The paper is organized as follows: In Section 2 we state the main results 
on convergence of the discretized Markovian iteration. The proof is given in 
several steps in Sections 3-5, where we establish the control of the Lipschitz 
constant, of the linear growth and the convergence of u^'"^ to respectively. 
In Section 6 we investigate the error due to the time discretization. To the 
best of our knowledge, our convergence theorem is the first of this type for 
coupled FBSDEs which also holds for a degenerate diffusion coefficient a. 
In Section 7 we spell out the proposed numerical scheme and present some 
numerical examples in Section 8. 



2. Notation and main results. The main results of this paper estimate 
the error of the discretized Markovian iteration (1.8) as the number of time 
steps n and the number of iterations m tend to infinity. Before we can state 
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these results, we need to fix some notation and discuss some assumptions. 
Prom now on we suppose, in the theoretical part, that all processes are one- 
dimensional. This is only to ease the notation and the attentive reader will 
easily see that all results hold true for the multi-dimensional case as well. 
The augmented filtration generated by the Brownian motion is denoted by 
F = {Tt,0<t<T}. 

The first assumption concerns the Lipschitz continuity and monotonicity 
of the coefficients. It will be in force throughout the whole paper without 
further notice. Denote 

Ax = xi- X2, Ay = yi-y2, Az = zi-Z2. 

Assumption 2.1. (i) There exist (possibly negative) constants kf,,kf 
such that 

[b{t,xi,y) - b{t,X2,y)]Ax < /c^lAxp, 
[f{t, X, yi,z)- f{t, X, 2/2, z)]Ay < /c/| Ay|^. 

(ii) b,a,f,g are uniformly Lipschitz continuous with respect to {x,y,z). 
In particular, there are constants K, by, ax, cry, fx, fz and gx such that 

\bit,xi,yi) - b{t,X2,y2)? < K\Ax\^ + by\Ay\^, 
\a{t,xi,yi) - fj(t,X2,y2)P < cTxlAxl"^ + aylAyl"^ , 
\f{t,xi,yi,zi) - f{t,X2,y2,Z2)\^ < MAx\^ + K\Ay\^ + f,\Az\\ 

lai^i) -g{x2)\^ < 5x|Axp. 

(iii) b(t, 0,0), a{t, 0,0), /(t, 0,0,0) are bounded. In particular, there are 
constants bo, a^, /o and go such that 

\b{t,x,y)\^ < bo + Klxl"^ + bylyl"^ , 

\a{t,x,y)\'^ < ao + crx\x\^ + dylyl"^ , 

\f{t,x,y, z)|2 <fo + fx\x\^ + K\y\^ + Mz\\ 

b(2;)l^ < 90 + 9x1x1"^ ■ 

We emphasize that here by et al. are constants, not partial derivatives. 
Indeed, we will not assume any differentiability conditions throughout this 
paper. For convenience, we also suppose that K is an upper bound for all 
the constants above. 

For results concerning the error due to the time discretization, we require 
the following assumption. 
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Assumption 2.2. The coefficients {b,a,f) are uniformly Holder-^ con- 
tinuous with respect to t. 

If Assumption 2.2 is in force, we use the same constant K to denote an 
upper bound of the square of the Holder constants. 

To ensure that a solution of (1.1) exists and the iteration converges, we 
further impose conditions which guarantee that we are in one of the following 
five cases: 

1. Small time duration, that is, T is small. 

2. Weak coupling of Y into the forward SDE, that is, by and ay are small. 
In particular, if by = ay = 0, then the forward equation in (1.1) does not 
depend on the backward one and, thus, (1.1) is decoupled. 

3. Weak coupling of X into the backward SDE, that is, fx and gx are small. 
In particular, if fx = 9x = 0, then the backward equation in (1.1) does 
not depend on the forward one and, thus, (1.1) is also decoupled. In fact, 
in this case Z = and (1.1) reduces to a decoupled system of a forward 
SDE and an ODE. 

4. / is strongly decreasing in y, that is, kj is very negative. 

5. b is strongly decreasing in x, that is, kb is very negative. 

The above conditions will be made precise later. 

Remark 2.1. We emphasize that Assumptions 2.1 and 2.2 alone are 
not sufficient to guarantee existence of a solution to the FBSDE (1.1). For 
instance, the one-dimensional linear FBSDE, 



does not admit a solution; see Ma and Yong [19] . 

We next present two examples from finance, in which we expect one of 
the conditions 1-5 to hold. For more details on both examples, we refer to 
Ma and Yong [19], Chapters 8.3 and 8.4. 

Example 2.1. (i) Consol rate models: In interest rate modeling the term 
structure can be determined by a so-called short rate r; see, for example, 
Brigo and Mercurio [5]. This short rate may depend on a long rate 
(also called consol rate), which in turn is influenced by the short rate via 
the relation 



dXt = Ytdt 



dVt 



Xt dt + Zt dWt 



t, 
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This problem can be cast into the FBSDE framework as 

drt = h{t, rt,Yt) dt + a{t, n^Yt) dWt, 

dYt = {rtYt-l)dt + ZtdWt, 

ro = R, Yt = 0. 

In generahzation of the Hull- White model [15] one can choose 

b{t, r, y) = K{e{t, y)-r), a{t, r, y) = i]{t, y) > 0. 

Here the drift h contains a mean reverting force k > 0, which pushes the short 
rate to a level which may depend on the long rate. The mean reverting force 
implies that b is monotonically decreasing in r, and so we expect condition 
5 to hold. 

(ii) Stock coupled with an option: In this example the forward SDE de- 
scribes a system of stock prices, and the backward SDE describes the price 
of an option on the stocks, leading in a complete market to the FBSDE 

dSt = bit, Yt)St dt + a{t, Yt)St dWt, 

dYt = {rYt + ZtOit, Yt)) dt + Zt dWt, 

So = s, YT = g{ST), 

where r denotes the riskless interest rate, 6 the premium of risk, and the 
drift b and volatility a of the stocks are influenced by the price Y of an 
option on S with pay-off function g. As the drift and the volatility of a 
stock fluctuate only slightly, we can expect that Y weakly couples into the 
forward part. 

Generically, we will derive the following theorems. The first theorem con- 
cerns the convergence of the iteration as m tends to infinity. 

Theorem 2.1. Under Assumption 2.1, let one of the conditions 1-5 
hold true. Then, for sufficiently small h, (1.6) has an ^^essentially" unique 
solution with linear growth and there are constants C > and < c < 1 
such that 

max Wl'^^ix) - < C(|x|2 + m)c^, 

0<i<n 

where u"'"* is given by (1.8). 

We will see from the proof that the constant c, which determines the 
rate of convergence, depends on the conditions 1-5. Roughly speaking, the 
stronger the monotonicity (resp. the weaker the coupling, the smaller the 
time horizon), the smaller one can choose c and, hence, the faster the itera- 
tion converges. 

Concerning the error due to the time discretization, we obtain the follow- 
ing: 
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Theorem 2.2. Suppose Assumptions 2.1, 2.2 hold true, and one of the 
conditions 1-5 is in force. Then equation (1-3) admits a viscosity solution 
u{t,x) with linear growth and there is a constant C > such that, for suffi- 
ciently small h, 



Note that the forward part in (1.6) is discretized by an Euler scheme. 
Hence, the stated convergence of order 1/2 for the time discretization is the 
best rate one can hope for. 

Combining these two theorems, one can derive the following with a little 
extra effort: 

Theorem 2.3. Under the assumptions of Theorem 2.2 FBSDE (1.1) 
has a unique solution [X, Y, Z) and there are constants C > and < c < 1 
such that, for sufficiently small h. 



These generic results will be made precise in Theorems 5.1, 6.3 and 6.5 be- 
low. We emphasize that none of the above theorems requires nondegeneracy 
of a and, in principle, X and W can have different dimensions. Moreover, 
we do not suppose any smoothness or boundedness conditions. However, we 
also underline again that FBSDE (1.1) does not allow coupling through the 
control part Z. 

The proof of convergence for the Markovian iteration, which will be given 
in Sections 3-5, is rather technical. We therefore briefly outline its proof. 

Strategy of proof. In a standard Picard iteration, like (1.7), one estimates 
|ym+i _ in terms of - X"'\ and then - X"'\ in terms of 

_ y™-i| . However, applying similar techniques to (1.8) yields only esti- 
mates of \X'^,rn+l _ xn,m^ ^^^.^^ |^n,m (j^n,m+l ) _ ^n,m-l ^^n,m^ | ^ gj^^^g 

yn,m _ ,m ^^^n ^ sccms Unavoidable to control the Lipschitz constant 
of n"'"' to obtain estimates in terms of |y"'™ — y">™-i|. 

More precisely, suppose, for the moment, that the step size n is fixed, 
and u"'"*(x) are bounded functions for all m and i. Then one can derive 



estmiates oi \X^ — | m terms oi sup^ \u^ [x) — {x)\ and a 



uniform (in time) Lipschitz constant L(u"''") of n"'™ by applying Lemma 



max |n"(x)-n(ti,x)p<C(l + |xp)/i. 

D<i<n 





<C(l + |x|2)[mc'" + /i]. 
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3.2 below. Combining this with estimates for sup^ \u^'"^~^^{x) — in 
terms of ^x^'^~^^ — X^'^\ (derived from Lemma 3.3 below), one obtains 

I n.m+l / \ n.m/ \\2 



(2.1) 



< c{L{u^n) sup |<''"(x) - <''"-'(x)p 



for some constant c(L(u"'"*)) which depends on the coefficients of the equa- 
tion and the Lipschitz constant of u"'™'. Now we wish to iterate the above 
estimate. To this end, we need a uniform control L of L{uP''^) and condi- 
tions on the coefficients which ensure that c(L) < 1. Section 3 is devoted to 
deriving such uniform control of the Lipschitz constants. 

In general, the conditions imposed in this paper do not guarantee that 
yn,m bounded, and therefore, the use of the sup-norm becomes meaning- 
less. However, one can easily see that n"'*" is of linear growth. Given linear 
growing functions there are constants G{(p) and H[Lp) such that 

W^{x)\'' <G{^)\x\'' + H{^) V(i,x). 

We consider linearly growing functions ip^ and Lp^ close to each other, if we 
can choose G{^p^ — 93^) and H{ip^ — Lp^) small. Following similar, but slightly 
more intricate considerations than the ones leading to (2.1), we can estimate 
Q{y,n,m+i _ yn,m) H{u'''"^+^ - n'^''") in terms of G(u"''" - u'^''""!) and 
H{u^'^ — ti"'*""^); see Theorem 5.2 below. However, the constant, which 
replaces the above c{L[vP'''^)) in these estimates, depends on the Lipschitz 
constant of n"'*" and additionally on the linear growth of tt"'"*~i through 
G(ii"'"*~^) and H{vP''^~^). Hence, for the general case a uniform control 
for the linear growth of tt"'™- is required as well. Such control will be given 
in Section 4. Then, iterating the above estimates yields the convergence 
of the Markovian iteration under each of the conditions 1-5, as will be 
demonstrated in Section 5. 

In order to study the behavior of the functions u"'"^, as outlined above, 
we introduce an important operator Fn for each n. For any measurable 
functions = {pi\Q<i<n-i-, define if) and $ as follows: 



(2.2) 



<^n{x)^g{x), 
[ Mx)^E{Y^4'^ + f{U,x,Y;^'i''',Mx))h} 



Xf^f = x + b{ti,x,pi{x))h + a{ti,x, pi{x))AWi, 



We finally set Fn{(p) = It is then obvious that u"'™ = F„(u"''^~^), and 
Fn{u^) = u"' if (1-6) has a solution u". We also point out that y"'™'^ given 
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by (1.8), can be expressed in the form 

J ti 

thanks to the martingale representation theorem. The analogous expression 
holds for y" defined in (1.4). 

3. Lipschitz continuity. In this section we obtain a Lipschitz constant of 
it"'™(x), uniformly in {i,n,m). To this end, we first investigate the Lipschitz 
continuity of Fnif). Given Lipschitz continuous ip, let L{(pi) denote the 
square of a Lipschitz constant of ipi, and L{ip) = supj -L((^j). Denote 

(3.1) 

Our aim is to derive the following theorem: 

Theorem 3.1. // 

(3.2) Lo<e-\ 

then for any L > Li and for h small enough, we have 

L(n"'™) < L Vm. 

Notice that (3.2) holds true in all five cases of Section 2. 
We prepare the proof of Theorem 3.1 with two lemmas. For constants 
Xj > 0, J = 1, 2, 3, denote 

Ai ^ 2kb + cj^. + 1 + Kh, 
A2 = by + ay + Kh, 

(3.3) ^3^A2 + A3 + (l + A2^)A'/i, 

^4 = 2kf + 1 + A3 V, + (1 + X2^)Kh, 
A5^f, + {l + \^^)Kh. 

Lemma 3.2. Fix i and for 1 = 1,2, let 

Xl, ^ X\ + h{ti, X\, ip'iXi))h + a{t„Xl,^\xl))AWi+u 

where X\ is !Fti-fn^O'Surable. Assume ip^ is uniformly Lipschitz continuous. 
Then for any Ai > 0, 

EtMXl, - Xl,\^} < [1 + A,h + (1 + Xi)A2hL{ip')]\Xl - Xf\^ 

+ {l + X^')A2h\ipHxf)-ip\xf)\'. 
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This lemma can be easily proved by some standard estimates and its proof 
is therefore omitted. 



Lemma 3.3. Fix i and for 1 = 1,2, let 



Y! = Yl^, + f(t„Xl,Yl^„ZDh- I ZldW^ 



where 



1 



Then for any A2, A3 > 0, 

I Ay,|2 + (1 - A3)/l|AZi|2 < (1 + A4/l)£iJ| Ay,+i|2} + ^5/l|AX, 

where 



2 



Proof. Denote 



^Z^Z'-Z\ Af ^ f{U,Xl,Y^_,„Zl) - fiU,Xf,Y^_,„Zf). 



Then 
(3.4) 



Ay, + / AZt dWt = AYi+i + Afh. 



Squaring both sides and taking conditional expectation, we have 



(3.5) 

Note that 



lAZA'^dt 



Et^{\AYi+i\' + 2AY,+iAfh + \Af\'h'}. 



Eu[£^' \AZt\^dtj > i Et,[£^' ^Ztdt 
By (3.4), we have 

Et^ I AZt dt\ = Et^ I AZt dWt AW,+i 

= EtA[^Y,+i + Afh]AWi+i} 
= h[AZi + EtAM^Wi+i}]. 
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Hence, 

> h\AZ^\^ + 2hAZ^Et^{AfAW^+l} 

> (1 - \2)h\AZi\^ - X2^h\Et^{AfAWi+i}\^ 

> (1 - A2)/i|AZ,|2 - X^^h^EtA\Af\^}. 
Thus, (3.5) implies that 

\AYi\^ + {1 - X2)h\AZi\^ 

(3.6) 

< Et,{\AY^+l\^ + 2AY+iAfh + (1 + A2 ')/i'| A/|2}. 
By Assumption 2.1(ii), we have 
(3.7) |A/|2 < K[\AX,\^ + |Ay,+i|2 + \AZ,W 

Moreover, 

where 



Af = A/i + A/2 + A/3, 



A/i — f{ti,Xl,Yj[^_i,Zl) - f{ti,X];,Yl^i,Zl), 

A/2 ^ /(t„ Xf, 1-1+1, Z/) - /(t,, XlYl^^, Zl), 

A/3 = /(ti,Xf ,yj^]^,Z/) - f{ti,Xf,Y■^_^_'^^,Zf). 
Then by Assumption 2.1(i) and (ii), we get 

2Al^;+iA/ = 2Ay,+iA/i + 2Ayi+iA/2 + 2Ay,+iA/3 
<|Ay,+i|2 + |A/i|2 + 2A:j|Ayi+i|2 

(3.8) + A3-V.|AK,+i|2 + A3/,-i|A/3p 

< |Ay,+i|2 + f,\AX,\^ + 2kf\AYi+i\^ 
+ A3V.|Ayi+i|2 + A3|AZi|2. 
Plugging (3.7) and (3.8) into (3.6), the lemma is proved. □ 

With these lemmas at hand, we can study the Lipschitz continuity of 
Fn{ip) given Lipschitz continuous ip. 

Theorem 3.4. For any Lipschitz continuous (p, we have 
L{Fn{ip)) < + A5T] [exp([^i + + AiA^hjT + [A2 + A2A4h]TL{p)) V 1] , 
where Ai = and A2, A3 > are chosen such that 

(3.9) ^3 < 1. 
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Proof. Recall (2.2). Fix i and xi,X2- Denote 

Ax = xi-X2, AX ^ x^'^'^i - Ay ^ y^'^'^i - 

^ ^>i(xi) - <i>iix2), A^Ji ^ iPiixi) - V'i(x2). 

We apply Lemmas 3.2 and 3.3, setting Ai = 0, and obtain 

E{\AXi+i\^} < [1 + Aih + A2hL{ip)]\Ax\'^ , 
|A$i|2 + (1 - ^3)/i| AV-iP < (1 + A4h)E{\AYi+i\'^} + As/ijAxp. 
By (3.9), we have 

|A$,|2 < [1 + Aih]L{<^i+i)E{\AX^+i\^} + A5h\Ax\^ 

< [1 + Aih][l + Aih + A2hL{ip)]L{^i+i)\Ax\'^ + ^/i] Ax|^. 

Thus, 

i(^0 < [1 + A^h] [1 + ^i/i + A2hL{yD)]L{<^i+i) + A^h 

(3.10) 

^ [1 + A/i]L($,+i) + As/i < [1 + A+h]L{^i+i) + A/i, 
where = ^ V and 

(3.11) A^Ai + Ai + AiAih + [A2 + A2Aih]L{^). 

Note that L(<&n) = S'x- Hence, we can apply the discrete Gronwall inequality 
to (3.10) and get 

Lm < e^^^ig. + A5T] = [g, + A,T][e'''' V 1], 
which, combined with (3.11), yields the assertion. □ 

We are now in position to give the proof of Theorem 3.1. 

Proof of Theorem 3.1. First, by induction, one can easily show that 
Lm — L(tt"'"*) < oo for each (n,m). Due to Theorem 3.4, we have 

Lm < [gx + A5T][eM[Ai +A4 + AiAih]T + [A2 + ^2^4/i]TL„_i) V 1], 

for Ai = and any A2, A3 > satisfying (3.9). 
Introducing 

L^^[A2 + A2Aih]TL^, 

we get 

Lrn < [^2 + A2^4/i]bx +^5r]r[e[^i+^4+^i^4'^]^e^™-i V 1] 

(3.12) 

< [A2 + A2Aih][g^ + ^5r]r[e[^i+^4+^i^4''l^e-^™-i + 1]. 
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Denote 

Lo{X,h) ^ [A2 + A2A4h][g, + A5T]Te^^'+^'^^''^^3.+A,^]TMAi+A,+Aa4h]T_ 
Obviously, Lq = 0. If 

(3.13) Lo(A,/i)<e-\ 
then, by induction, one can easily show that 

Lm<[A2 + A2A4h][g^ + A5T]T+l Vm. 
We plug this into the right-hand side of (3.12) to obtain 

Lm < [A2 + A2A^h][g, + ^5r]r[e[^i+^*+^i^4''l^+[^2+^2^^'^l [^-+^5^1^+1 V 1]. 
Thus, 

Lm < [9x + A5r][e[^i+^4+AiA4fc]T+[^2+A2A4h][9.+A5T]T+l ^ 

(3.14) 

= LiiX,h). 

So we want to choose A2, A3 and h which satisfy (3.9) and minimize Lq{X, h). 
Recall again that Ai = 0. In dependence of h we set, for small /i, 

(3.15) \2{h)^Vh, X^{h)^l-[l + K]Vh-Kh. 

Then ^3 = 1 and 

limLo(A(/i), /i) = Lq, limLi(A(/i), h) = Li. 

hlO hlO 

Suppose that (3.2) holds true. Then for any L > Li, we obtain Lo(A(/i), h) < 
and Li{X{h),h) < Z, provided h is small enough. In view of (3.14), the 
theorem is proved. □ 

4. Linear growth. This section is devoted to studying the linear growth 
of the functions ^"'"^(x). Given linearly growing functions ipi, assume 

\ip^ix)\^<G{ip^)\x\^ + H{ip^) VX, 

and let 

G{lp) = sup G{Lpi), H{(p) = sup H{(pi). 

i i 

To state the main result of this section, we first introduce the functions 

(4.1) ro(x)^^^, T^{x,y)^ sup ee'^Toidy) Vx,yGM; 
X o<e<i 
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and for G > 0, 

co(G) ^ T[g,Ti{[2kf + 1 + f,]T, [2kb + 1 + a,]T + [by + ay]GT) 

+ /.rro([2A:/ + 1 + /,]r)ro([2A;b + 1 + a,]T + [by + ay]GT)]; 
ci{G)^[by + ay]co{Gy, 

L2{G) ^ et^'^z+i+Z^l+^^o + foTro{[2kf + 1 + f,]T) + [6o + aoMG). 

Theorem 4.1. Assume (3.2) holds true and 

(4.2) ci(Li)<l. 

For any G > Li, ci(Li) < ci < 1, L2 > L2(Li), and for h small enough, we 
have 

G('u"''") < G, H(u''''^) < Vm. 

1 - Cl 

Notice that 

lim ro(x) = 0, 

x— »— 00 

(4.3) lim Ti{x,y) = 0, 

X—f — OO 

lim ri(x, y) = 0. 

Hence, (4.2) is satisfied in cases 1-5 of Section 2. 

Again we start with some a-priori estimates whose proofs are fairly straight- 
forward and hence omitted. Denote 

^1 — ^0 + o"o + Kboh, 

(4.4) 

B2^fo + Kfoh. 

Lemma 4.2. Assume 

Xi+i =Xi + b{ti,Xi,ip{Xi))h + a{ti,Xi,ip{Xi))AWi+i. 

Then, 

EtA\X,+i\^} < [1 + Alh + A2hG{ip)]\Xi\^ + [Si + A2H{^)]h. 
Lemma 4.3. Assume 

Yi = Yi+i + f{ti, Xi,Y,+i,Zi)h - Zt dWt, 

Jti 

where 

Zi = ]-Eu{Yi+lAW^+l}. 
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Then, for any A2, A3 > 0, 

|y,|2 + (1 - A3)h\Z,\^ < [1 + Aih]Et^{\Y^+^\^} + A^hlXil^ + B2h. 

To derive bounds for the linear growth of Fn{<f), we define discrete time 
versions of Fq and Ti by 



^ {1 + xhy-i 



(4.5) 

TUx,y)^ sup (l + xhmiy) 

0<i<n 

and discrete time versions of co{G),ci{G), L2{G) by 

co(A, h, G) ^ 5xr? (A4, Ai + ^2G) + A5r^(^4)r^(^i + a^g), 

(4.6) ci(A,/i,G)^^2Co(A,/i,G), 

L2(A, /i, G) ^ i?iCo(A, /i, G) + [e^4^ V l]go + S2rS(^4). 



Theorem 4.4. For ant/ linearly growing ip, 

(4.7) G(F„(99)) < [5,. + ^5r][e[^i+^^+^i^4'^]^+[^2+^2^4'^l^^(^) V 1], 

(4.8) H{Fnm < ci(A, /i, G((/.))/7(v9) + L2(A, h, G{^)), 
where A2, A3 > are supposed to fulfill (3.9). 

Proof. Denote ^ = Fn{^). Fix (zq, x) and define, for i = i^, . . . ,n — 1, 



(4.9) 



^Xi + h{ti,Xi,^i{Xi))h + a{ti ,Xi,^i {Xi ) ) A Wi+i , 
Yn^g{Xn), 

Yi = Y,+i + f{ti, Xi,Yi+i,Zi)h - [ 

Jti 



ZtdWt. 



Obviously Yi^ = <I>j„(a;). We obtain from Lemma 4.2 that 

E{\X^+i\^} < [1 + Aih + A2/iG((/p)]F{|X,|2} + [Si + A2Hi^)]h. 

Then 

i?{|X,|2} < [1 + Aih + A2/iG((^)]*-^oi5;{|X,„ |2} 

+ [Bi + A2H{ip)]h ^ [1 + A^h + A2hG{^)]^-'° 

(4.10) 



li-«0 U|2 



[l+^i/i + ^2/iG(v?)]*-*°|x| 
+ + A2H{^)]Ti~''^ {Ai + ^2G((^)). 
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Next, applying Lemma 4.3 and by (3.9), we have 

E{\Yi\''} < [1 + Aih]E{\Yi+^\''} + A^hE{\X,\''} + ^3/1. 

Note that 

\Yn? < go+9x\Xn\'^- 

Then 

\'^ioix)\^ = \Yio\^ < {l + A^hr-'<>[go + gxE{\Xn\^}] 



+ A,h ^ (1 + A^hy-'<^E{\Xi\^} + B2T'^-'"{A^). 



I=l0 



This, together with (4.10), imphes 
G{^J < (1 + A^hr~''^g,[l + Aih + A2hG{^)r-'<^ 

+ E (1 + A^hy-''' [1 + Aih + A2hG{ip)Y-'° , 



l=«0 

H{^i,) < (1 + ^4/1)"-^" 50 + i32r^'»(^4) 



+ [Bi + A2H{cp)] 



5.(1 + A4/i)"-^«r^-^«(Ai + A2G{^)) 

n-l 

+ A,hJ2i^ + A^hy-'"Ti-'° (Ai + A2Gi^)) 



Note that, for <i <n, 

(l + x/i)*<e^^Vl, rUx) <n(x), 

(4.11) 

{l + xhyTliy)<TUx,y). 

Then 

G($,„) < [^^, + ^5r][e[^i+^^+^i^^'^]^+[^2+^2^4'^l^^(^) V 1], 
H{^i,) < [e^*^ V l]5o + 52r^(^4) + co(A, h, G)[Bi + A2i?((^)]. 
Since the right-hand side does not depend on io, the assertion is proved. □ 

After these preparations we give the proof of Theorem 4.1. 

Proof of Theorem 4.1. Denote Gm = G{u'^'"'),Hm = i?(u"'™). Ob- 
viously, Go = Hq = 0. We may now conclude from Theorem 4.4 that, under 
(3.9), 

(4.12) Gm < [gx + A5r][e[^i+^4+^i^*''l'^+[^2+A2A4h]TG„-i y 

(4.13) Hrn < Cl(A, h, Gm-l)Hm-l + L2{X, h, Gm-l). 
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We now choose X-sih) and X4{h) as in (3.15) for small h. Since (3.2) holds 
true, for any G > Li, we can follow the same arguments as in Theorem 3.1 
and get G(u"'™) < G from (4.12). Note that 

jim^r^(x) =rro(xr), jim r^(x,y) = rri(xr,yr), 

limci(A(/i),/i,G) =ci(G), limL2(A(/i), /i, G) = L2(G). 

hlO hlO 

For any ci, ci(Li) < ci < 1, and L2, L2{Li) < L2, we can choose G > Li 
such that ci(G) < ci and L2{G) < L2. Then, for sufficiently small h, it holds 
that ci{\{h),h,G) <ci and L2{\{h),h,G) < Lg. Now, by (4.13), we get 

Hm < CiHm-1 + L2, 

which implies the result. □ 

5. Convergence of the Markovian iteration. We now make the assump- 
tions of Theorem 2.1 precise and prove convergence of the Markovian itera- 
tion as the number of iteration steps tends to infinity. 

To this end, we first introduce 

C2(Ai, L, G) = [el^k,+i+a^+[by+ay]G]T V 1] (1 + X~^)[by + ay]T 
X [g,Ti{[2kf + l + f,]T, 

[2h + 1 + CT^. + (1 + Ai)[6j, + ay]L]T) 
+ f^Tro{[2kf + l + f,]T) 

X ro([2A;b + 1 + a:, + {1 + Xi)[by + ay]L]T)]; 
C2{L,G)= mi C2(Ai,L,G). 

Ai >0 

We will prove the following theorem: 

Theorem 5.1. Assume (3.2) and 
(5.1) C2(Li,Li)<l. 



(i) For any L > Li,G > Li, L2 > L2{Li),ci{Li) < ci < 1, there exists a 
solution n" to (1.6) such that 

(5.2) L{u'')<L, G{u'')<G, H{u'^)<H^ 



1-ci' 
if h is small enough. 
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(ii) For any C2{Li,Li) < C2 < 1 and for h small enough, 

I n,m/ \ n/ \\2 

max |u. • {x)-u,{x)\ 



(5.3) 




3 



— + [bo + ao + {by + ay)H]TG mc^. 



+ 



(1 - v^2r 



(iii) Fix G > and suppose is another solution to (1.6) with linear 
growth such that G{vP') < G. Then vP' = u^, if h (depending on G) is small 
enough. 

Remark 5.1. (i) In view of (4.3), it is straightforward to see that (5.1) 
is also satisfied in cases 1-5 of Section 2. 

(ii) One can check directly that ci(L) < C2(L, L), and thus, condition (5.1) 
implies (4.2). 

Again we first study the operator F„ to prepare the proof of Theorem 5.1. 

Theorem 5.2. Assume V"'^)^^ have linear growth and (p^ is Lipschitz 
continuous. Then, for any Ai > 0, 



G(F„(^i)-F„(v.2)) 

<C2{Xi,h,L{ip'),G{ip^))G{ip' - ^^), 
H{Fn{^^)-Fn{v^)) 

<C2{Xi,h,L{ip'),G{ip^))H{ip' -ip^) 

+ C2{\i,h, L{^^),G{ip^))[Bi + A2H{ip^)\TG{^^ - y.^), 



Proof. For Z = 1, 2, denote = Fn{ip^). Fix (zq, x) and define {X\Y\ P) 
analogously to (4.9). Then obviously y/^ = ^{^{x). Denote L = L{ipi), and 
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Application of Lemma 3.2 yields 

£;{|AXi+i|2} 

(5.5) 

< ^{[1 + A^h + (1 + Xi)A2hL]\AXi\^ + (1 + X:[^)A2h\A^{Xf)\^}. 

Note that 

\A^{Xf)\' <G{A^)\Xf\' + H{A^). 
By the first inequality in (4.10) and (4.11), 

(5.6) sup E{|X2|2}<[|x|2+[Si + A2i/(<^2)]2^][g[Ai+A2G(v2)]T^-^] 

io<i<n 

Then (5.5) implies 

E{\AXi+i\^} < [1+Aih + (1 + Xi)A2hL]E{\AX^\^} 
+ (1 + X^^)A2h[G{A^)A + H{Aip)]. 

Since AXi^ = 0, we get 
sup ^{ I AX, 1 2} 

io<i<n 

< (1 + Ar^)^2/i[G'(A(^)i + H{Aip)] 

(5.7) 

n— 1 



^[l+^l/l+(l + Ai)A2/ii]*"^° 



X 



= (1 + Xi^)A2[G{A^)A + H{A^)]T^-'%Ai + (1 + Ai)^2L). 
Furthermore, we obtain from Lemma 3.3 and (3.9), 

E{\AY,f}<[l + A4h]E{\AYi+i\^} + Ar,hE{\AXi\^}. 
Hence, by (5.7) and (4.11), 

mMf = \AY,f 

<(l + A4/i)"-^°i^{|Ay„|2} + ^5r^'-'o(A4) sup E{\AX,\^} 

iQ<i<n 

<[(l + ^4/i)"-*°5x + A5^^'"(^4)] sup E{\AXi\^} 

iQ<i<n 

< (1 + X^')A2[g.{l+A,hr-'<^r'S-'%Ai + (1 + Xi)A2L) 

+ (^4)r^*" (^1 + (1 + Xi)A2L)] 

X [G(A(^)i + /^(Ay?)] 

< (1 + X^^)A2[G{Aip)A + F(A(/j)] 
X [g,TUA4,Ai + {l + Xi)A2L) 

+ A5T^{A4)T^{Ai + (1 + Ai)^2i^)], 



MARKOVIAN ITERATION FOR COUPLED FBSDES 21 
which, together with (5.6), imphes the theorem. □ 

We can apply this theorem to estimate the distance between u"'"* and 

yn,m—l 

_ Theorem 5.3. Assume that L(u"''") < L, G(u"''") < G and < 
H for all m£N and sufficiently small h. Moreover, assume 

C2{L,G)<1. 

Then for any C2{L,G) < C2 < 1 and for h small enough, 

(5.8) ^(n^'^+i - n"'™) < Gc^; 

(5.9) /7(n"'™+i - n"'"^) < + [(6o + <ro) + {by + ay)H]TGm]c^. 

Proof. Choose A2,A3 depending on h as in (3.15). Note that with this 
choice 

hm r^{x, y) = TTi{xT, yT), hm C2(Ai, /i, L, G) = C2(Ai, L, G). 

/i— »0 n— >0 
Hence, we may find an appropriate Ai such that, for h small enough, 

C2{Xi,h,L,G) < C2, 
C2{\i,h,L,G)[Bi + A2H] < C2[[bo + c7o] + [by + ay]H]. 

Denote 

Applying Theorem 5.2, we get, for small h, 

(5.10) G(An"''"+i) < C2G(An'^''"), 

(5.11) i/(Au"'™+i) < C2if(Au"'™) + C2[6o + CTO + [by + ay]H]TG{Au''''^). 
Note that 

G(An"'i) = G(n"'i)<G, 

By (5.10), we get (5.8). Moreover, together with (5.8), (5.11) implies (5.9) 
immediately. The proof is complete now. □ 

Theorem 5.1 can now be proved by iterating the above theorem. 

Proof of Theorem 5.1. Assume G, L, L2, ci, C2 satisfy the conditions 
specified in the theorem. Without loss of generality, we assume C2{L, G) < 02- 
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Recall that (5.1) implies (4.2). Hence, by Theorems 3.1 and 4.1, we get, for 
h small enough, 

L(n"'™) < L, G(n"'™) < G, //(u"'"*) < H. 

Hence, (i) will follow directly from (ii). 
To prove (ii), we denote 

L^[[bo + ao] + [by + ay]H]TG. 

Applying Theorem 5.3, we get 

K'^'^Hx) - <'™(x)p < [Gjxp + H + Lm]c^. 

Then 

|<''"+^(x) - <''"(x)| < [^/Q\x\ + Vn+ 4im]^^'\ 
Thus, for any mi > m. 



''"(x)-nr^(x)|<5: + 



m 



1- V m (1-^)2 

Note that the right-hand side above converges to as m — > oo. Then u^'^{x) 
is a Cauchy sequence and hence converges to some u^{x). Moreover, 



K''"(x) - n7(2;)|2 < 3[G\x\^ + H + Lm] 



.m 



which leads to (5.3) and thus proves (ii). 

It remains to prove (iii). For any G > 0, assume is another solution to 
(1.6) with hnear growth such that G{u'^) < G. Then Fn{u^) = n". Note that 
tt" = = n". Assume u^^-^ = u^-^-i- We now apply a local version of Theorem 
5.2. That is, we consider (2.2) only on the interval with terminal 

condition $i+i(x) = ttf^i(x) (instead of on [0,T] with terminal condition 
g{x)). We note that in this case there is only one time subinterval. One can 
check directly that 

Tl{x) = h, T\{x,y) = {l+xh)h. 

Then (5.4) becomes 

C2{h) ^ [e[^i+^2G]h V 1](1 + X-^)A2[L{1 + A4h)h + A.h"]. 
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Set if^ = u^,ip'^ = u". We note that, for given ip, G{ip) and H{ip) are not 
unique. So, in general. Theorem 5.2 does not lead to 

G«-nr)<C2WG(nr-tir), 

- u^i) < C2{h)H{vJi - + C2{h)[Bi + A2H{u'')]TG{u'l - n^). 

To get around this difficulty, we use the definitions of G{ip), H{(p). Let Gq,Hq 
be some constants satisfying 

72^(^)1' <G'o|x|2 + Fo. 

For = 1, 2, . . . , denote 

Now applying Theorem 5.2 repeatedly for z/ = 1, 2, . . . , we get 
(5.12) K(x)-n^(x)p<G^|x|2 + i7^ 

Note that 

= Go52(/^)^ = HQC2{hY + [Si + ^2ii'(u")]Tz/C2(/i)^ 

For h small enough, we have C2{h) < 1. Then 

lim Gu = lini H,y = 0. 

Sending — > oo in (5.12), we get n" = uf. Repeating the arguments back- 
war dly, we prove that = u^. □ 

6. Convergence of the time discretization. We now study the error due 
to the time discretization. We first introduce a continuous time version of 
the operator Fn- Suppose 99 is a function on [0, T] x R which is Lipschitz in 
the space variable and let (X'^''"'^', y^'*^''^, Z'''''^''^) be the unique solution to 
the decoupled FBSDE (0 < r < t < T) 

' ^V'-nx = x + j' b{s, X,^'^'^ y?(s, Xl'''-'^)) ds 

+ jy{s,Xf'^'-,^{s,Xf'^n)dWs, 

y^'^'"'^ = gix^'"'"") + r f{s, y/''^'^ ''■'^) ds - ^ z^^"" dw^. 
^ Jt Jt 

We then define ^{t,x) = Y-^'^'^ and F{{p) = It is known from Pardoux 
and Peng [21] that, under Assumption 2.2 and if (p is additionally continuous 
as a function in time and space, $ is a viscosity solution to the following 
semilinear PDF: 

r $f + io-2(t, X, (f)^^^ + b{t, X, (f)<^^ + f{t, X, $^a(t, x,ip)) = 0, 
U(T,x) = <7(x). 
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We now define recursively vP = and n™ = Then the fohowing 

theorem can be proved similarly to, actually more easily than, Theorem 5.1. 
A detailed proof can be found in the appendix of the preprint version, which 
is available from the authors upon request. 

Theorem 6.1. Assume (3.2) and (5.1) hold true. 

(i) u"* converges to some function u uniformly on compacts. 

(ii) \u{t,xi) - n(t,X2)|2 < Lilxi - xap; \u{t,x)\'^ < Li|x|2 + j^l^. 

(iii) \u{t, x) — u{s, x)p < C(l + |xp)|t — s\ for some constant C . 

(iv) F{u) = u. Moreover, if F{u) = u and u has linear growth, then u = u. 

(v) Under Assumption 2.2, u is a viscosity solution to (1.3). 

From now on we denote by C a generic constant which may depend on 
the coefficients b,a,f,g, but is independent of n,/i and x. The value of C 
may vary from line to line. 

We next consider the following decoupled FBSDE: 

rt rt 



(6.1) 



Xt = x+ f b{s,Xs,u{s,Xs))ds+ f a{s,Xs,u{s,Xs))dWs 
Yt = g{XT)+j^ fis,Xs,Ys,Zs)ds- Z^dWs, 



and its time discretization: 
' X^ = X, 



(6.2) <^ 



Xr+i4xr + 6(t,,Xr,n(ti,Xr))/i + a(ti,Xr,n(ti,Xr))At^,+i, 



-n I ' 

rn A 



Denote nO(x) = u{t^,x) and = F„(u°). It is obvious that = <(Xf ). 
Note again that (6.1) is decoupled. By Theorem 6.1, and applying nowadays 
standard arguments for decoupled FBSDEs (see, e.g., Delarue [9] and Zhang 
[26] for (i), and Zhang [25] and Bouchard and Touzi [4] for (ii)), we can derive 
the following corollary. A detailed proof is again given in the appendix of 
the preprint version. 

Corollary 6.2. Assume all the conditions in Theorem 6.1 hold true. 

(i) FBSDE (1.1) has a unique solution {X,Y,Z), which also solves (6.1), 
and it holds that Yt = u{t,Xt). 

(ii) Moreover, we have the following estimates: 

(6.3) |n7(x) - u{ti,x)\^ < C{1 + |x|2)/i, 
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sup e\ sup [|Xt-lf„i|2 + |yt-y^'i,|2]\ 

l<i<n Ue[t,_i,t,] J 

(6.4) +x:i^{/j^t-^r-iprfi] 



i=l 

<C7(l + |x|2)/i. 

Applying the above decoupling relation, Yj = u{t,Xt), and the conver- 
gence results for decoupled FBSDEs, stated in (6.3)-(6.4), we can establish 
the convergence of u", as the time grid becomes finer. 

Theorem 6.3. Suppose Assumption 2.2 is in force, and (3.2) and (5.1) 
hold true. Then 

lu'lix) - u{ti,x)\^ < C[l + \x\^]h. 

Proof. For any L > ci(Li) and G > ci(Li), when h is small, we have 
L{u'')<L, G{u'')<G. 
Moreover, we know from Theorem 6.1(ii) that 

L(n°)<Li, G('u°)<Li, H{u°)<H ^ ^ 



l-ci(Li) 

Note that Fn{u^) = u". Applying Theorem 5.2 on u" and , we get 
G('«" - u") < C2(Ai,/i,Z,G')G(u" - -u"), 
//(u*" - w") < C2(Ai, /i, L, G)F(u" - n°) 

+ c2(Ai, /i,Z, G)[Si + ^2^]rG(n" - n°). 
For any e > 0, we obtain, thanks to (6.3), 

|<(x) - n0(x)|2 < (1 + e)|<(x) - u^{x)\^ + C,|nr(x) - n(t„x)|2 

< (1 + e)[G(n" - {i")|x|2 + i/(n" - -S")] + Q(l + |x|2)/i 

< [(1 + e)c2(Ai, /i, Z, G)G(n" - n°) + Ce/i]|x|2 
+ (l + e)c2(Ai,/i,Z,G) 

X [i7(n" - u°) + [Si + A2H]TG{u'' - u^)] + Ceh. 

Now for any 02(^1, -Li) < C2 < 1, we can choose Z,^ and e appropriately 
such that, for h small enough, 

(1 +e)c2(Ai,/i,Z,G) < C2. 
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Then we get 

- n°(x)|2 < [c2G(tx" - n°) + Ceh]\x\^ 

(6.5) 

+ C2//(n" - n°) + C7eG(n" - n°) + C7^/i. 

We now follow the arguments in the proof of Theorem 5.1(iii). Fix some 
Go, Hq such that 

K(x)-7xO(x)|2<Go|x|2 + //o. 

For = 1, 2, . . . , denote 

Gy = C2Gy-l + Ceh, Hy = C2Hy-l + G^Gy-l + G . 

Then (6.5) implies that 

(6.6) K(x)-n°(x)|2<G^|x|2 + i?^ "iu. 

Note that 



1-C2' 



1-C2 



1-C2 



1/C2 



1 -c^ 
+ G,/i- ^ 

1 -C2 



Since C2 < 1, and sending ^ 00 in (6.6), we get 

iur(x)-u°(x)|2<-^N2 + -^. 

1 - C2 (1 - C2Y 

The proof is complete. □ 

As a direct consequence of Theorems 5.1 and 6.3, we have the following: 

Theorem 6.4. Under the assumptions of Theorem 6.3 we have, for any 
€2(^1, Li) < C2 < 1 and for h small enough, 

K'™(x) - u(ti,x)p < C(l + |x|2)[mc^ + h]. 

We close the theoretical part of this paper with a precise version of the 
generic Theorem 2.3. 

Theorem 6.5. Under the assumptions of Theorem 6.3 we have, for any 
C2{Li,Li) < C2 < 1 and for h small enough, 

sup e{ sup [|Xi-xf_:7|2 + |y,-y^:^'™|2]\ 

l<i<n Ue[ti-i,ti] J 

+ E^{/' \Zt-Zj'l'^\^dt^<Gil + \x\^)[mc^ + h]. 
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Proof. By (6.4), it suffices to prove 

n-1 

(6.7) sup E{\AXi\'^ + \AYi\'^} + hY,E{\AZi\'^}<C{l + \x\'^)[mc^ + h], 

0<i<n 

where 

First, by Lemma 3.2 with Ai = 1, we get 

E{\AXi+if} < E{{1 + Ch)\AXif + C/i|n(ti,Xf •'") - <''"(Xf ''")|2} 
< (1 + Ch)E{\AXi\^} + C{1 + \x\^)[mdp + h]h. 
Since AXq = 0, we have 

J2 E{\AXi\^}<C{l + \x\^)[mc^' + h]. 

Q<i<n 

Next, choose A2 = A3 = and h smah enough so that A3 < ^. Applying 
Lemma 3.3, we obtain 

^{lAy^l^ + i/i|AZ,|2} < ^{(1 + Ch)\AYi+i\^ + Ch\AX^\^}. 

Since 

|Ay„|2 = \giX^) - 5(X-'-)|2 < C\AXn\\ 

we can easily get 

n-1 

sup E{\AY,\^} + hJ2E{\AZi\^} 

0<i<n 

<C sup E{\AXi\^}<C{l + \x\'^)[mdp + h]. 

0<i<n 

This proves (6.7) and hence the theorem. □ 

7. A numerical algorithm. We now briefly explain how the discretized 
Markovian iteration above can be transformed into a numerical algorithm 
which is viable also for high-dimensional problems. To this end, we replace 
the conditional expectations by a simulation based least squares regression 
estimator, as was suggested, for example, by Gobet, Lemor and Warin [14] 
and Bender and Denk [3] in the context of decoupled FBSDEs. An alter- 
native estimator based on Malliavin calculus is discussed in Bouchard and 
Touzi [4] for decoupled FBSDEs. A quantization algorithm for reflected BS- 
DEs is presented in Bally and Pages [2]. 

For the reader's convenience, we spell out our algorithm for the coupled 
case. While a convergence analysis is out of the scope of the present paper. 
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we will illustrate the algorithm by some numerical examples in the next 
section. 

We assume that the number of time steps n is fixed for the remainder of 
this section. In the algorithm conditional expectations are first replaced by 
orthogonal projections on K basis functions. Then the orthogonal projec- 
tions are approximated by simulating A trajectories. Hence, the algorithm 
can be described for the one-dimensional case iteratively as follows. It is 
straightforward how this extends to the multi-dimensional case: 

• Fix some xq. Set ul'^'^'^{x) = 0. 

• Sample A independent copies of the time discretized Brownian motion 
W^:, i = 0, . . . ,n, A = 1,...,A, starting in and denote the corresponding 
increments by AW^. 

n -n,m—l,K,A/ \ ■ i i . j. j t j. v-n,m.\ A i 

• Suppose (x) IS already constructed. Let Aq = xq and 

where — for notational convenience — we suppress the dependence of x"'"^'^ 
on K through ^^"."^-1'^.^. Note, depends on all Brownian incre- 

ments AW^, i = 1, . . . ,n, A = 1, . . . , A, through ^j^'™"^'-^'^. While we ex- 
pect that this dependence will make a convergence analysis difficult, the 
examples below indicate that the algorithm works without re-simulating 
the Brownian paths in every iteration step. 

• Choose a set of Lipschitz continuous basis functions 

such that 

(7.1) ^^n,m,k^j^n,rn,X^^ k=l,...,K} 

forms a subset of L^(r2). Prom the construction below, it will become 
evident that u^'"^'^'^{x) inherits the Lipschitz continuity from the basis 
functions. This feature seems to be important to ensure that the dis- 
cretized forward equations for x"'™^^'^ do not explode. 

• Define, for i = n— 

-n,m,K,Af ^\ A „/\ ;r.n,m,K,A/ ^\ A 




v^' ' ' (x) = argmf<j -y^+'i ' AW^+i - ' ) 



2 
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V G span(i3f ''"'^) 



I A=l 



Note that the minimization problems are hnear least squares problems, 

which can be easily implemented. 

Let 

Yi =ui' ' (Xi' • ), 

1^1 - - 

A=l 



1 ^ 

-0-71,171, K,X A ^ v>n,m,ii',A , -r>n,m,K,X F^n,m,K,X\j 

^0 =X^-^i +/(0,xo,ri' ,^0 

A=l 



We expect that the thus constructed (X"'"*'^, y".™.-?^.^, ^n,m,i^,A) ^^^.^ "^lose" 
to (X"'™'-^, y"'.m'.^^^"."i.->^)^ the solution of the discretized Markovian it- 
eration (1.8) with the Brownian motion W replaced by W^, if the basis 
functions are chosen appropriately and the number A of simulated paths is 
sufficiently large. While an analysis of the error by estimating the condi- 
tional expectations is left to future research, the numerical examples in the 
next section support this conjecture. 

8. Numerical examples. For the simulations, we consider the example 
Xd,t = Xd,o + [ cfYu dWd,u, 

3 

Yt = Y,^iTv{Xd,T)+ r-ry„+ie-3^(^-"V2('^sin(X,,„)) du 
d=i \d=i / 

J* d=i 

where Wd,t^ d=l,...,D, isa D-dimensional Brownian motion and a > 0, r, 
Xd,o are constants. Note that the corresponding differential operator degen- 
erates at y = 0. 
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By Ito's formula, one can easily check that this FBSDE decouples via the 
relation 

D 

(8.1) yt = e-^(^-*)5]sin(Xrf,i). 

d=l 

Note that for small a the weak coupling condition of Y into X is satisfied, 
while, for large a, the monotonicity condition of / can be fulfilled by choosing 
r large enough. 

In the simulations we replace conditional expectations by least squares 
regression as explained above with the "canonical" basis functions 

1, Xd, l<d<D, {-R)V (xdXq) AR, l<d<q<D, 

that is, monomials up to order two in x = (xi, . . . ,X£)). The truncation con- 
stant R guarantees that the basis functions are Lipschitz continuous. We 
set 

R^W, Xd,o = ^, l<d<D, 

r = l, A = 50000, n = 50, 

unless otherwise stated. With this initial condition, we get Yq = De~^'^-^~^^ . 
Recall also that the estimator y^"'™-'-'^'-^ Qf Yq does not depend on A and is 
denoted by Yq'"^'^ from now on. 

Figure 1 illustrates the convergence of the iteration in the case of a four- 
dimensional state space (Z? = 4). Both figures display the absolute error 
|y^n,m,x _ ^ function of the number of iterations m. In Figure 1(a) 

the case r = (no monotonicity) is considered for several values of a which 
represent different influences of the coupling. In Figure 1(b) the coupling 
parameter a = OA is fixed, while the strength of the monotonicity varies 
by different values of r. In general, we observe that the iteration converges 
extremely fast, as could be expected in view of Theorem 2.1 which states 
mc*", for some c, < c < 1, as rate of convergence. From the proof of this 
theorem we know that c is the smaller, the weaker the coupling or the 
stronger the monotonicity is. This explains the faster convergence observed 
for small values of a and large values of r. 

The influence of the time partition is displayed in Figure 2. It shows 
the absolute error |y^"'"^='°P'^ — ^o| as a function of the number of time 
points n. We stop the iteration when two consecutive estimates Yq'"^'^ are 
within a distance of lO"'^. This iteration level is denoted mstop- The observed 
convergence rate is in accordance with l/\/n as derived in Theorem 2.2. 

Finally, we demonstrate that the space dimension four is no limitation 
for the proposed algorithm. To this end, we consider the lO-dimensional 
case with the parameter values r = 0, o" = 0.1 in Figure 3. Under the same 
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Absolute error as lunction of the number of iterations 




number of iterations 

(b) 

Fig. 1. Convergence of the iteration for different choices of cr and r. 
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Absolute error as the number of time steps increases 

0.1 6 I 1 r 1 1 1 1 1 r 




number oi lime steps 

Fig. 2. Convergence as the time partition becomes finer for r = 0, a = 0.4, D = A. 



stopping criterion as above, the iteration terminates after 12 steps. Recall 
that, from (8.1), we can also approximate the true value of X via the usual 
Euler scheme (applying the same simulated Brownian increments l^W-^). 
The corresponding approximation along the Ath path is denoted X"'"^, and 
hence, 

D 



d=l 

may be considered a close approximation of Yf.. In Figure 3 we display, 
for the 10-dimensional case, a comparison between a typical path of 1^"'^° 
(dashed line) and Y:"-'™-st°P'^''*'o (^^^{[(^ line), as well as the (absolute) empirical 
mean square error between Y^''^ and y^"''"^*"?' ' ^ A = 1,...,A. Precisely, 
Figure 3(b) shows 



1 A 

A l-'i ^i \ 

A=l 
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as function of time. 

Remark 8.1. Figure 3 shows a larger mean square error close to ter- 
minal time than close to initial time. This is a rather typical feature, when 
conditional expectations are estimated by the above least squares method. It 
can be explained by interpreting this method as creating a stochastic mesh 
(see Glasserman [13]) and observing that this mesh is typically much finer 
close to initial time than close to terminal time. Hence, the error close to ter- 
minal time, observed in Figure 3, cannot be diminished by solely increasing 
the number of iterations, but by improving the quality of the conditional ex- 
pectation estimator. A generic trick to improve the quality close to terminal 
time is to add the terminal function g to the basis. 

Acknowledgment. The paper benefited from the constructive remarks of 
an anonymous referee. 
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